REPORT  DOCUMENTATION  PAGE 


form  Approved 
0MB  No.  07Q4-01B8 


aath«finS»ndmjintatnlfw  th«  daU  nei^ed.  and  completing  and  reviewing  the  collection  of  information.  Send  comments  reoardlng  this  bor^  estimate  v  any  thh 

informatloninduding  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  Information  Operatiiw  and  RefKKts,  1215  iefferson 
S2^H%hwi?.  sISS^  the  Office  of  Managenfent  and  Budget.  Paperwork  Reduction  Prelect  (0704^J18«).  Washington.  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  b!ank)  2.  REPORT  DATE 

1997 


4.  TITLE  AND  SUBTITLE 

Matrix  methods  for  population  analysis 


3.  REPORT  TYPE  AND  DATES  COVERED 

TECHNICAL 


5.  FUNDING  NUMBERS 


NSF  DEb  92-11945 
NSF  OCE  93-02874 
ONR  N00014-92-J-1527 


7.  PERFORMING  ORGANIZATION  NAME(S}  AND  AOORESS(ES) 


WOODS  HOLE  OCEANOGRAPHIC  INSHTUTION 
WOODS  HOLE,  MA  02543 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


WHOI  CONTR.  9224 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADORESS(ES)  10.  SPONSOWN^MONITORING 

AGENCY  REPORT  NUMBER 

OFFICE  OF  NAVAL  RESEARCH 
ENVIRONMENTAL  SCIENCES 
DIRECTORATE 

ARLINGTON,  VA  22217-5660 _  _  _ 


11.  SUPPLEMENTARY  NOTES  ,  ,  -  i.  ij  t. 

In  citing  this  report  in  a^bibliography,  the  reference  given  snould  be. 

In:  S.  Tuljapurker  and  H.  Caswell.  1997.  Structured-population  models  in  marine,  terrestrial 

and  freshwater  systems.  Chapman  and  Hall,  New  York. 


12a.  DISTRIBUTION /AVAILABIUTY  STATEMENT 

APPROVED  FOR  PUBUC  RELEASE: 
DISTRIBUTION  UNLIMITED 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 


iwnx 


14.  SUBJECT  TERMS 


population  model s 
age  structure 
size  structure 


15.  NUMBER  OP  PAGES 

40 


16.  PRICE  CODE 


17  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  I  19.  SECURITY  CLASSIFICATION  I  20.  LIMITATION  OF  ABSTRACT  ■ 

OF  REPORT  i  OF  THIS  PAGE  I  OF  ABSTRACT  I  I 

UNCLASSIFIED  1  UNCLASSIFIED  |  UNCLASSIFIED  | 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 


S.  Tuljapurkar  and  H.  Caswell.  1997. 
Structured-Population  Models  in  Marine, 
Terrestrial  and  Freshwater  Systems. 
Chapman  and  Hall,  New  York. 


CHAPTER  2 

Matrix  Methods  for 
Population  Analysis 


Hal  Caswell 


Matrix  models  for  structured  populations  were  introduced  by  P.  H. 
Leslie  in  the  1940's  (Leslie  1945,  1948).  Although  they  are  in  some 
ways  the  simplest  of  the  mathematical  approaches  to  structured 
population  modeling  (see  Chapter  1),  their  analysis  requires  com¬ 
putational  power.  For  this  reason,  and  because  ecologists  of  the  day 
viewed  matrix  algebra  as  an  esoteric  branch  of  advanced  mathe¬ 
matics,  they  were  largely  neglected  until  the  late  1960’s,  when  they 
were  rediscovered  by  ecologists  (Lefkovitch  1965)  and  human  de¬ 
mographers  (Goodman  1967;  Keyfitz  1967).  In  the  1970’s,  matrix 
models  were  adopted  by  plant  ecologists,  who  discovered  that  they 
could  easily  handle  the  complexity  of  plant  life  cycles  in  which  size 
or  developmental  stage  was  more  important  than  chronological  age 
in  determining  the  fate  of  individuals  (Sarukhan  &  Gadgil  1974; 
Hartshorn  1975;  Werner  &  Caswell  1977). 

This  chapter  introduces  the  construction  and  analysis  of  matrix 
population  models.  I  will  not  try  to  be  comprehensive;  I  have  done 
that  elsewhere  in  book  form  (Caswell  1989a)  and  twice  in  simplified 
form  with  a  focus  on  particular  taxa  (Caswell  1986;  McDonald 
&  Caswell  1993).  Instead,  I  try  to  convey  the  basics  of  matrix 
population  models  clearly  and  briefly.  Wherever  possible,  I  use 
different  derivations  than  before  (Caswell  1989a),  so  you  may  find 
some  new  ways  to  understand  the  source  of  some  familiar  results. 
I  rarely  cite  my  book  (Caswell  1989a)  (in  spite  of  having  done 
so  three  times  in  this  paragraph);  almost  every  topic  presented 
here  could  be  followed  by  the  instruction,  “see  the  book  for  more 
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information.”  My  focus  here  is  on  methods;  I  am  sparing  in  my 
use  of  examples,  because  they  can  be  found  in  many  of  the  other 
chapters  in  this  volume. 

A  note  about  notation.  I  use  boldface  symbols  to  denote  vectors 
(lower  case,  as  in  n)  and  matrices  (uppercase,  as  in  A).  Entries 
of  vectors  and  matrices  are  lowercase  letters  with  subscripts,  so 
the  zth  entry  of  n  is  n^,  and  the  element  in  the  fth  row  and  jth 
column  of  A  is  a^j.  Sometimes  I  use  parenthetical  superscripts  to 
label  matrices  or  vectors.  Thus,  Am  or  A^”^^  might  both  be  used 
to  denote  the  mth  in  a  series  of  matrices;  the  ijth  element  of  this 
matrix  is  written  The  transpose  of  the  matrix  A  is  A^.  If 
X  =  a  +  is  a  complex  number,  the  complex  conjugate  is  denoted 
hy  X  =  a  —  bi.  The  complex-conjugate  transpose  of  A  is  A*.  The 
scalar  product  of  two  vectors  is  (x,  y)  =  y*x. 


1  Formulating  Matrix  Models 

A  matrix  population  model  operates  in  discrete  time,  projecting  a 
population  from  t  to  t  +  1.  The  first  step  in  formulating  a  matrix 
model  is  to  define  the  time  scale  for  the  projection;  this  is  called  the 
projection  interval  Models  for  the  same  population  with  different 
projection  intervals  may  look  quite  different. 

The  second  step  is  to  choose  a  set  of  state  variables  for  individ¬ 
uals  {i-state  variables):  these  provide  the  information  necessary  to 
determine  the  response  of  an  individual  to  the  environment,  over  a 
projection  interval.  Examples  of  i-state  variables  include  age,  size, 
developmental  stage,  and  geographical  location. 

A  matrix  model  uses  discrete  stages,  so  the  third  step  is  to  define 
a  set  of  discrete  categories  for  each  f-state  variable.  Some  i-states 
are  naturally  discrete  (e.g.,  instars),  while  others  are  naturally  con¬ 
tinuous  and  must  be  made  discrete  (e.g.,  size).  Dividing  continuous 
variables  into  discrete  categories  involves  trade-offs.  A  model  treats 
all  individuals  within  a  category  as  identical,  so  creating  only  a  few 
large  categories  reduces  the  accuracy  of  the  z-state  dynamics.  Cre¬ 
ating  many  small  categories,  alternatively,  leads  to  a  large  model 
and  may  make  it  hard  to  estimate  parameter  values  because  sample 
sizes  in  each  category  are  small. 

The  stages  describe  the  life  cycle,  or  as  much  of  it  as  we  believe 
to  be  demographically  important.  The  next  step  is  to  translate 
them  into  a  model.  The  life-cycle  graph  is  a  useful  tool  for  this 
translation. 
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The  Life- Cycle  Graph 

A  life-cycle  graph  describes  the  transitions  an  individual  can  make, 
during  a  projection  interval,  among  the  ^-state  categories  that  de¬ 
fine  its  life  cycle.  To  construct  the  graph,  first  draw  a  numbered 
point  (a  “node”  in  graph- theory  terminology)  for  each  i-state  cat¬ 
egory.  If,  for  example,  size  is  the  z-state  variable,  then  the  life-cycle 
graph  contains  a  node  for  each  size  class.  If  age  and  size  are  both 
i-state  variables,  then  the  life-cycle  graph  contains  a  node  for  each 
age-size  category.  Draw  arrows,  “directed  arcs,”  between  nodes  to 
indicate  where  it  is  possible  for  an  individual  in  one  stage  to  con¬ 
tribute  individuals  to  another  stage  over  a  single  projection  inter¬ 
val.  The  head  of  the  arrow  shows  the  direction  in  which  individuals 
move.  If  individuals  can  contribute  in  both  directions  between  two 
stages,  draw  two  arrows,  rather  than  an  arrow  with  a  head  on  both 
ends.  Contributions  from  one  stage  to  another  can  result  from  the 
movement  of  individuals  from  one  stage  to  another  (e.g.,  by  growth 
or  aging)  or  from  production  of  new  individuals  (e.g.,  by  birth). 

With  each  arrow  is  associated  a  coefficient;  the  coefficient  on 
the  arrow  from  stage  j  to  stage  i  is  denoted  aij  (the  ordering  of 
the  subscripts  is  important;  it  corresponds  to  the  arrangement  of 
coefficients  in  the  resulting  matrix  model).  The  coefficient  aij  gives 
the  number  of  stage  i  individuals  at  t  + 1  per  stage  j  individual  at 
time  t. 

So  far,  we  have  made  no  decisions  about  the  nature  of  these 
coefficients;  I  return  to  this  below. 

Figure  la  shows  a  life-cycle  graph  for  an  age-classified  model 
with  the  age  interval  equal  to  the  projection  interval.  Individuals 
in  one  age  class  can  contribute  to  another  only  by  surviving  to  the 
next  older  age  class  or  by  reproduction  to  the  first  age  class.  Figure 
16  shows  the  graph  for  a  size-classified  model  in  which  individuals 
may  grow  to  the  next  size  class,  remain  in  their  own  size  class,  and 
possibly  reproduce  new  individuals  in  the  first  size  class.  Suppose 
some  individuals  in  the  first  size  class  grow  so  rapidly  that  after 
one  projection  interval  they  are  in  the  third  size  class.  This  would 
require  the  modification  shown  in  Figure  c,  1  where  an  arrow  has 
been  drawn  from  stage  1  to  stage  3. 

The  interpretation  of  the  coefficients  depends  on  the  identity  of 
the  stages  and  the  processes  involved  in  the  transitions.  In  Fig¬ 
ure  la,  the  Pi  are  age-specific  survival  probabilities  and  the  Fi  are 
age-specific  fertilities.  In  Figures  16  and  Ic,  the  Gi  are  probabili¬ 
ties  of  surviving  and  growing,  the  Pi  are  probabilities  of  surviving 
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Figure  1.  Three  life-cycle  graphs,  (a)  An  age-classifed  model  with  three 
age  classes;  the  Pi  are  age-specific  survival  probabilities,  and  the  Fi  are 
age- specific  fertilities,  (b)  A  size- classified  model  with  three  size  classes; 
the  Gi  are  size-specific  probabilities  of  survival  and  growth,  the  Pi  are 
size-specific  probabilities  of  surviving  and  remaining  in  the  same  size 
class,  and  the  Fi  are  size-specific  fertilities,  (c)  The  same  life-cycle  graph 
as  in  (b);  but  with  an  additional  transition  (Hi)  from  size  class  1  to  size 
class  3. 

and  not  growing,  and  the  Fi  are  size-specific  fertilities.  In  Figure 
Ic,  the  coefficient  Hi  is  the  probability  of  surviving  and  growing 
enough  to  move  from  size  class  1  to  size  class  3,  Demographers 
use  the  term  vital  rates  to  refer  collectively  to  the  rates  of  sur¬ 
vival,  growth,  reproduction,  and  any  other  important  demographic 
processes. 


A  Set  of  Difference  Equations 

The  life-cycle  graph  corresponds  directly  to  a  model  written  as  a  set 
of  difference  equations.  For  the  size-classified  graph  in  Figure  lb, 
remembering  the  definitions  of  the  coefficients,  the  set  of  equations 
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describing  the  population  is 

ni[t  +  1)  =  Pini{t)  +  F2n2{t)  +  Fznz{t) , 

+  1)  =  Gini(t)  +  P2n2{t) ,  (1) 

n^[t  +  1)  =  G2n2{t)  +  Pznz{t) . 

It  is  worth  looking  at  these  equations  for  a  moment.  The  first  states 
that  the  number  of  individuals  in  stage  1  at  t  +  1  is  the  sum  of 
those  remaining  in  stage  1  from  time  t  and  those  contributed  by 
reproduction  from  stages  2  and  3.  The  second  equation  states  that 
the  number  in  stage  2  at  t  +  1  is  the  sum  of  those  growing  into 
stage  2  from  stage  1  and  those  remaining  in  stage  2  from  time  t. 
The  third  equation  says  the  same  thing  for  stage  3. 

The  equations  corresponding  to  Figure  Ic  are 

ni(t  +  1)  =  Pini{t)  +  F2n2{t)  +  Fznz{t) , 

n2(t-fl)=  P2n2{t)  ^  (2) 

Tizit  +  1)  =  H-  G2'^2{^)  'b  • 

It  would  be  possible  to  write  down  these  equations  directly,  with¬ 
out  using  the  life-cycle  graph,  if  we  were  clear  about  the  nature  of 
the  possible  transitions,  which  in  turn  depends  on  the  definition  of 
the  stages.  But  using  the  life-cycle  graph  makes  it  easier,  and  helps 
to  guard  against  mistakes  in  defining  the  stages  and  transitions. 


The  Matrix  Model 


The  system  of  difference  equations  derived  from  the  life-cycle  graph 
can  be  written  more  simply  in  matrix  form: 


n{t  +  1)  =  An(t) , 


(1) 


where 


n(t)  = 


ni(t) 

n2{t) 


nk{t) 


(2) 


is  a  stage- distribution  vector  and  A  is  a  population-projection  ma¬ 
trix.  The  elements  of  this  matrix  can  be  obtained  from  the  system 
of  difference  equations  or  directly  from  the  life-cycle  graph:  the 
ijth  entry  of  A  is  the  coefficient  on  the  arrow  from  stage  j  to 
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stage  i.  The  reason  for  the  order  of  the  subscripts  is  to  guarantee 
this  correspondence. 

Applying  this  rule  to  the  life-cycle  graphs  in  Figure  1  yields 
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The  age-classified  model  produces  a  special  matrix,  with  positive 
entries  only  on  the  first  row  (fertilities)  and  the  subdiagonal  (sur¬ 
vival  probabilities).  Such  a  matrix  is  often  called  a  Leslie  matrix, 
in  recognition  of  the  early  papers  of  Leslie  (1945,  1948). 

I  have  said  nothing  about  how  the  numerical  values  of  the  coeffi¬ 
cients  aij  are  obtained.  This  obviously  important  question  deserves 
its  own  chapter  (see  Chapter  19,  by  Wood,  for  one  approach),  but 
here  I  assume  that  the  matrix  is  at  hand  and  focus  on  how  to 
analyze  it. 


Types  of  Matrix  Models 

The  coefficient  Oij  is  the  contribution  of  each  individual  in  stage 
j  to  the  number  of  individuals  in  stage  i  during  one  projection 
interval.  What  happens  in  the  next  projection  interval?  Depending 
on  the  answer  to  this  question,  matrix  models  fall  into  three  classs, 
each  with  its  own  analytical  approach. 

Linear,  constant- coefficient  models.  If  the  coefficients  aij  are  con¬ 
stants,  the  resulting  model  is  linear  and  time-invariant: 

n{t  -h  1)  =  An(t) .  (6) 

This  is  the  simplest  case:  it  can  be  analyzed  in  great  detail,  and 
it  is  widely  used.  But  in  reality  the  vital  rates  are  not  constants, 
so  the  biological  interpretation  of  these  results  requires  great  care. 
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NouliueaT  models.  If  the  aij  are  not  constant  but  depend  on  the 
current  state  of  the  population,  the  resulting  model  is  nonlinear: 

n(t  +  1)  =  Ann(t) ,  (7) 

where  An  is  the  transition  matrix  evaluated  at  n.  The  nonlin- 
earity  may  result  from  density  dependence  (e.g.,  competition  for 
resources),  frequency  dependence  (e.g.,  competition  for  mates), 
or  both. 

Time-varying  models.  The  coefficients  may  also  change  with  time, 
indcpondoiitly  of  n(f).  resulting  in  the  model 

n(^  +  1)  =  Atn{t).  (8) 

l)rTrTininisTi(\  periodic  variation  is  often  used  to  model  seasonal- 
iT\-  or  otlKT  kinds  of  environmental  periodicity.  Alternatively,  the 
coefficients  may  vary  stochastically,  reflecting  some  random  envi¬ 
ronmental  process  (see  Chapter  3,  by  Tuljapurkar).  Time- varying 
models  may  be  either  linear  or  nonlinear. 

Objectives  of  Analysis 

The  analysis  of  each  of  these  types  of  model,  although  requiring 
different  mathematical  tools,  addresses  a  set  of  similar  questions. 
Imagine  that  you  are  in  possession  of  a  matrix  population  model. 
What  you  should  do  with  it  depends  on  the  question  you  want  to 
answer. 

1.  Transient  analyses  describe  the  short-term  dynamics  resulting 
from  specific  initial  conditions. 

2.  Asymptotic  analyses  describe  the  long-term  dynamics  of  the 
population. 

(a)  Population  growth  rate:  what  is  the  asymptotic  rate  of  pop¬ 
ulation  growth  or  decline? 

(b)  Population  structure:  what  are  the  relative  abundances  of  the 
different  stages  in  the  life  cycle? 

(c)  Ergodicity:  are  the  dynamics,  including  the  growth  rate  and 
the  population  structure,  asymptotically  independent  of  ini¬ 
tial  conditions? 

(d)  Attractors  (mainly  in  density-dependent  models):  what  are 
the  qualitative  properties  of  the  asymptotic  dynamics  (fixed 
point,  cycle,  quasiperiodicity,  chaos,  etc.)? 
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3.  Perturbation  analyses  examine  the  effects  of  changes  in  param¬ 
eter  values  or  initial  conditions  on  the  results  of  the  models. 
Three  questions  are  of  particular  importance: 

(a)  Sensitivity  and  elasticity  analysis  of  population  growth  rate: 
how  does  the  growth  rate  respond  to  changes  in  vital  rates? 

(b)  Stability  analysis  of  equilibria:  if  initial  conditions  are  per¬ 
turbed  slightly  away  from  an  equihbrium  point,  does  the 
solution  return  to  or  depart  from  the  neighborhood  of  the 
equilibrium? 

(c)  Bifurcation  analysis:  what  happens  to  the  asymptotic  be¬ 
havior  of  a  nonlinear  model  as  a  parameter  in  the  model  is 
changed? 

The  methods  used  to  address  these  questions  depend  on  the  nature 
and  sometimes  on  the  details  of  the  model,  but  any  population¬ 
modeling  project  that  does  not  address  short-term  dynamics,  long¬ 
term  dynamics,  and  the  effects  of  perturbations  on  those  dynamics 
has  left  something  out. 

2  Analysis:  The  Linear  Case 

We  begin  with  the  linear  time-invariant  model  (6),  in  which  A  is 
a  constant  matrix.  There  are  two  justifications  for  spending  time 
on  this  model,  in  spite  of  the  fact  that  the  vital  rates  of  any  real 
population  are  certainly  not  constant.  The  first  is  theoretical:  un¬ 
derstanding  population  dynamics  in  the  simplest  case  is  a  funda¬ 
mental  step  in  understanding  more-complicated  cases.  The  second 
is  practical:  when  interpreted  as  a  projection  rather  than  a  predic¬ 
tion  (Ke5ditz  1968;  Caswell  1989a),  the  results  of  a  linear  model 
provide  a  valuable  characterization  of  the  current  environment  by 
calculating  the  purely  hypothetical  consequences  of  maintaining 
that  environment  forever.  Linear  matrix  population  models  are 
frequently  used  in  this  way,  as  a  form  of  demographic  analysis 
of  vital-rate  data,  rather  than  as  a  prediction  of  future  population 
dynamics. 


Exponential  Solutions  and  the  Characteristic  Equation 

One  approach  to  equation  (6)  is  to  conjecture  that,  like  other  linear 
equations,  it  has  an  exponential  solution(s), 

n(t)  =  A^w 


(9) 
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for  some  fixed  vector  w.  Substituting  this  into  (6)  gives 
\t+i^  rr  A^Aw. 

A  scalar  A  and  a  vector  w  that  satisfy  this  relation  are  called  an 
eigenvalue  and  eigenvector  of  A,  respectively.  They  must  satisfy 

(A  -  AI) w  =  0 , 

which  has  a  nonzero  solution  for  w  only  if  the  determinant  of  the 
matrix  A-AI  equals  zero.  This  is  called  the  characteristic  equation: 

det(A-- AI)  =0.  (10) 


The  Spectral  Decomposition  of  A 

An  alternative  approach  is  to  begin  by  solving  equation  (6),  start¬ 
ing  from  a  specified  initial  population  n(to)-  By  repeatedly  apply¬ 
ing  (6),  we  see  that  n{to  4- 1)  =  An(to)?  n(^o  +  2)  =  A^n(to),  3.nd 
in  general 

n(to  + 1)  =  A*n(to)  •  (11) 

Thus,  to  understand  population  dynamics  over  time  we  need  only 
understand  the  behavior  of  A*. 

One  approach  to  the  problem  is  via  the  spectral  decomposition  of 
A,  which  makes  it  possible  to  evaluate  any  function  of  A,  including 
A^  First,  note  a  few  facts  about  the  eigenvalues  and  eigenvectors 
of  a  matrix.  The  vectors  w  and  v  are  right  and  left  eigenvectors  of 
A  if  there  is  a  (possibly  complex)  scalar  A  such  that 

Aw  =  Aw,  (12) 

v*A  =  Av*,  (13) 

where  the  asterisk  denotes  the  complex-conjugate  transpose.  A  left 
eigenvector  v  of  A,  corresponding  to  A,  is  a  right  eigenvector  of  A* 
corresponding  to  A;  that  is, 

A*v  =  Av.  (14) 

The  eigenvalues  are  found  as  the  solutions  of  the  characteristic 

equation  (10). 

If  A  is  a  fc  X  fc  matrix,  the  characteristic  equation  is  a  polynomial 
of  degree  fc  and  has  fc  solutions  A^,  i  =  1, 2, . . . ,  fc.  The  correspond¬ 
ing  eigenvectors  are  w^  and  v^,  i  =  1,2,. I  assume  that  these 
eigenvalues  are  all  distinct,  as  seems  to  be  true  in  practice  for 
population-projection  matrices.  This  assumption  guarantees  that 
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the  right  eigenvectors  and  left  eigenvectors,  respectively,  are  lin¬ 
early  independent  sets. 

Let  (w,  v)  =  v’" w  denote  the  scalar  product  of  w  and  v.  The  left 
and  right  eigenvectors  can  always  be  scaled  so  that  (v^,  w^)  =  1.  In 
addition,  the  left  and  right  eigenvectors  corresponding  to  different 
eigenvalues  are  orthogonal,  so  that  (vj,  Wj)  =  0  if  z  j. 

Any  matrix  A  with  distinct  eigenvalues  can  be  written  in  the 
form 


A  —  AiZi  H - 1-  AfcZjt , 

where  the  matrices  Z^,  known  as  the  constituent  matrices  of  A,  are 
given  by 

Zi  =  .  (15) 

That  is,  Zi  is  a  matrix  whose  columns  are  all  proportional  to 
and  whose  rows  are  all  proportional  to  v*. 

The  constituent  matrices  have  two  important  properties.  First, 

Zf  =  Wivjw^v* 

=  Wi(wi,Vi}v*  (16) 

-  Z,. 

(Such  matrices  are  called  idempotent.)  Second,  multiplying  two 
different  constituent  matrices  yields  a  zero  matrix: 

ZiZj  =  ^iV-WjVj 

=  Wi(w^-,Vi)v;  (17) 

=  0. 

These  properties  are  useful  because,  together,  they  imply  that 

A2  =  AjZj'j  =  ^  A?Z, .  (18) 

Multiplying  repeatedly  by  A,  it  is  not  hard  to  see  that 

A*  =  ^A%.  (19) 

1 

This  result,  together  with  equation  (13),  yields  our  desired  expres¬ 
sion  for  the  dynamics  of  a  population  described  by  (8); 

n(to  +  t)  =  ^AfZin(to). 

2 


(20) 
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The  only  parts  of  the  right-hand  side  of  (20)  that  vary  with  time 
are  the  factors  Af .  The  behavior  of  Af  depends  on  the  sign  of  Xi  and 
on  whether  Aj  is  real  or  complex.  If  Xi  is  real  and  positive,  Aj  grows 
or  decays  exponentially,  depending  on  whether  Xi  is  greater  or  less 
than  one.  If  Xi  is  real  and  negative,  A|  oscillates  between  positive 
and  negative  values,  growing  or  deca}dng  in  magnitude  depending 
on  whether  |Ai|  is  greater  or  less  than  one.  If  Aj  is  complex,  A| 
oscillates  in  a  sinusoidal  pattern,  growing  or  decaying  in  magnitude 
depending  on  whether  |Ai|  is  greater  or  less  than  one. 

Since  the  dynamic  properties  of  the  population  are  determined 
by  the  eigenvalues  of  A,  it  behooves  us  to  see  what  we  can  say,  a 
priori  and  in  general,  about  these  eigenvalues. 


Eigenvalues,  Eigenvectors,  and  the  Perron-Frobenius  Theorem 

We  can  safely  assume  that  the  elements  of  A  are  nonnegative.  Neg¬ 
ative  elements  in  A  imply  the  possibility  of  negative  individuals, 
which  I  prefer  not  to  deal  with.  Perhaps  surprisingly,  this  sim¬ 
ple  assumption  tells  us  almost  everything  we  want  to  know  about 
the  eigenvalues  and  eigenvectors  of  A,  thanks  to  a  mathematical 
result  known  as  the  Perron-Frobenius  theorem.  In  order  to  state 
the  theorem,  we  need  two  more  properties  of  A:  irreducibility  and 
primitivity. 

A  matrix  A  is  irreducible  if  and  only  if  its  life-cycle  graph  is 
connected,  that  is,  if  there  is  a  path,  following  the  direction  of 
the  arrows,  from  every  stage  to  every  other  stage.  A  matrix  A  is 
primitive  if  and  Only  if  there  is  some  integer  k  such  that  every 
element  of  A*=  is  strictly  greater  than  zero.  A  more  biologically 
revealing  criterion  is  based  on  the  life-cycle  graph.  Define  a  loop  as 
a  sequence  of  arrows,  traversed  in  the  direction  of  the  arrows,  that 
begins  and  ends  at  the  same  node,  without  passing  through  any 
node  twice.  The  matrix  A  is  primitive  if  and  only  if  the  greatest 
common  divisor  of  the  lengths  of  the  loops  in  the  life-cycle  graph 
is  one.  Any  primitive  matrix  is  also  irreducible.  Most  population- 
projection  matrices  encountered  in  practice  are  both  irreducible 
and  primitive. 

What  about  matrices  that  are  reducible  or  imprimitive  (i.e.,  not 
primitive)?  A  reducible  matrix  has  some  stages  that  make  no  con¬ 
tribution  to  some  other  stages;  the  life-cycle  graph  breaks  into 
two  (or  more)  pieces  with  only  one-way  communication.  The  most 
common  example  is  a  life  cycle  with  post-reproductive  stages;  from 
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Figure  2.  (a)  A  life-cycle  graph  corresponding  to  a  reducible  matrix. 
Stage  4  represents  post-reproductive  females;  there  is  no  pathway  from 
this  stage  to  any  of  the  earlier  stages,  (b)  A  life-cycle  graph  correspond¬ 
ing  to  an  imprimitive  matrix.  This  is  a  semelparous  age-classified  model; 
individuals  that  survive  to  age  class  3  reproduce  and  die. 

such  a  stage  there  is  no  pathway  back  to  the  part  of  the  life  cycle 
that  does  reproduce.  Figure  2a  shows  an  example;  a  graph  of  this 
form  appears  in  a  stage-classified  model  for  killer  whales  (Brault  & 
Caswell  1993).  An  imprimitive  life  cycle  has  some  underlying  pe¬ 
riodicity,  so  that  the  loops  in  the  life-cycle  graph  are  all  multiples 
of  some  common  loop  length.  Imprimitive  matrices  are  sometimes 
called  “cyclic”  to  reflect  this  fact.  The  most  common  example  is  a 
semelparous  age-classified  life  cycle  with  a  fixed  age  at  reproduc¬ 
tion  (Fig.  26).  Only  a  single  loop  appears  in  such  a  life-cycle  graph, 
with  a  length  determined  by  the  age  at  reproduction.  Some  kinds 
of  seasonal  models  for  annual  organisms  also  produce  imprimitive 
matrices,  reflecting  the  periodicity  imposed  by  the  annual  cycle  of 
the  seasons.  The  graphs  in  Figure  2  contain  no  coefficients  because 
reducibility  and  primitivity  depend  on  the  form  of  the  graph  but 
not  on  the  values  of  the  coefficients. 

The  Perron-Probenius  theorem  states  that  a  nonnegative,  irre¬ 
ducible,  primitive  matrix  has  three  properties: 

1.  a  simple  (i.e.,  non-repeated)  eigenvalue  \i  that  is  real,  positive, 
and  strictly  greater  in  magnitude  than  any  of  its  other  eigenval¬ 
ues, 

2.  a  right  eigenvector  Wi  corresponding  to  Ai,  which  is  strictly 
positive  (or  can  be  made  so  by  multiplying  by  a  scalar)  and  is 
the  only  nonnegative  right  eigenvector,  and 

3.  a  left  eigenvector  Vi  corresponding  to  Ai,  which  is  also  strictly 
positive  and  is  the  only  nonnegative  left  eigenvector. 
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The  Perron-Probenius  theorem  also  describes  the  eigenvalues  and 
eigenvectors  of  imprimitive  and  reducible  matrices  (see  Caswell 
1989a  and,  of  course,  many  matrix  texts,  e.g.,  Seneta  1981;  Horn 
&  Johnson  1985). 

Demographic  Ergodicity 

The  eigenvalue  Ai  (the  dominant  eigenvalue)  plays  a  central  role  in 
the  asymptotic  analysis  of  linear  matrix  models.  In  equation  (20), 
the  growth  of  n{to  + 1)  is  given  by  a  sum  of  terms  involving  the 
eigenvalues  of  A  raised  to  higher  and  higher  powers.  Intuitively,  as 
t  gets  large,  Aj  increases  more  quickly,  or  decreases  more  slowly, 
than  Af  for  i  ^  1.  Asymptotically,  we  expect  the  growth  of  the 
population  to  be  determined  by  Aj,  whereas  all  the  eigenvalues 
contribute  to  short-term  transient  behavior.  More  precisely, 

limy^f-^^  Zin(to) 

Zin(to)  (21) 

wiVin(to)  • 

This  gives  the  following  results  on  asymptotic  dynamics,  condi¬ 
tional  on  the  primitivity  of  A, 

1.  The  population  eventually  grows  geometrically  at  a  rate  given 
by  Ai  (the  population  growth  rate  or  rate  of  mcreose). 

2.  Population  structure  eventually  becomes  proportional  to  wi 
(the  stable  stage  distribution). 

3.  The  constant  of  proportionality  relating  population  structure 
and  w  is  a  weighted  sum  of  the  initial  numbers  in  each  stage 
(v];n(to))*  The  weights  are  the  elements  in  vj;  the  vector  vi 
thus  gives  the  relative  contributions  of  the  stages  to  eventual 
population  size  {not  population  growth  rate)  and  is  called  the 
reproductive-value  vector. 

The  population  eventually  converges  to  the  stable  stage  distribu¬ 
tion,  growing  at  a  rate  given  by  the  dominant  eigenvalue,  regardless 
of  the  initial  conditions  (except,  of  course,  the  special  case  of  a  zero 
initial  population) .  The  property  of  forgetting  the  past  and  grow¬ 
ing  at  a  rate  determined  by  the  vital  rates  rather  than  by  initial 
conditions  is  called  ergodicity. 

Because  Ai,  wi,  and  vj  are  properties  of  the  vital  rates  rather 
than  initial  conditions,  they  are  widely  used  as  demographic  statis- 
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tics.  They  can  provide  valuable  insight  into  the  vital  rates  and  the 
environmental  conditions  that  determine  them,  but  they  cannot 
predict  actual  population  dynamics;  everyone  “knows”  that  popu¬ 
lations  cannot  grow  geometrically  forever.  These  statistics  must  be 
interpreted  as  projections  of  what  would  happen  if  the  vital  rates 
were  to  remain  constant,  rather  than  as  predictions  of  what  will 
happen.  They  characterize  the  present  environment,  not  the  future 
of  the  population. 

Similar  ergodic  results  hold  for  stochastic  models  (see  Chapter 
3)  and  density-dependent  models.  In  each  case,  the  asymptotic 
properties  provide  demographic  statistics  that  are  determined  by 
the  vital  rates  rather  than  by  the  historical  accidents  of  initial 
conditions.  They  can  be  used  just  as  Ai,  Wi,  and  Vi  are  used  in 
the  linear  case. 


3  Perturbation  Analysis 

Only  rarely  are  we  interested  in  one  precisely  specified  model.  We 
can  usually  imagine  that  the  model  might  change  in  some  way, 
and  would  like  to  know  how  such  changes  would  affect  the  re¬ 
sults  of  the  analysis.  Perturbation  analyses  address  this  problem. 
In  density-independent  models,  perturbation  analyses  focus  on  the 
eigenvalues  and  eigenvectors,  whereas  in  density-dependent  mod¬ 
els  perturbation  analyses  focus  on  the  stability  and  bifurcation  of 
equilibria. 

A  perturbation  analysis  of  the  eigenvalues  of  a  population-pro¬ 
jection  matrix  can  answer  several  questions. 

1.  What  are  the  effects  of  potential  changes  in  the  vital  rates, 
as  might  result  from  strategies  designed  to  protect  endangered 
species  (by  increasing  A)  or  control  pest  species  (by  reducing 
A)? 

2.  Where  should  efforts  to  improve  the  estimates  of  the  vital  rates 
be  focused  in  order  to  improve  the  accuracy  of  the  estimate  of 
A?  All  else  being  equal,  the  biggest  payoff  comes  from  improving 
the  estimates  of  the  vital  rates  to  which  A  is  most  sensitive,  since 
errors  in  those  estimates  have  the  biggest  effect. 

3.  Genetic  variation  produces  individuals  whose  vital  rates  are  per¬ 
turbed  from  the  overall  population  values;  from  these,  natu¬ 
ral  selection  chooses  those  perturbations  whose  carriers  increase 
most  rapidly.  Which  vital  rates  are  imder  the  greatest  selective 
pressure? 
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4,  Suppose  that  some  environmental  differences  (either  natural  or 
the  result  of  experimental  manipulation)  have  produced  differ¬ 
ences  in  the  vital  rates,  and  hence  in  A,  among  two  or  more 
populations.  How  much  do  each  of  the  vital-rate  differences  con¬ 
tribute  to  these  observed  differences  in  A? 

Fortunately,  it  is  easy  to  calculate  the  sensitivity  of  A  to  a  change 
in  any  of  the  vital  rates,  once  we  know  the  eigenvectors-  The  next 
subsection  presents  these  calculations.  Formulas  also  exist  for  the 
sensitivities  of  the  eigenvectors  w  and  v,  for  the  sensitivities  of  A 
for  periodic  time- varying  models  (Caswell  Trevisan  1994),  and 
for  the  sensitivities  of  the  sensitivities  themselves  (Caswell  19966). 


Sensitivity  and  Elasticity  of  Eigenvalues 

The  sensitivity  of  population  growth  rate  to  changes  in  the  vital 
rates  can  be  calculated  as  the  derivative  of  A  to  changes  in  the 
matrix  elements  aij.  Suppose  that  A,  w,  and  v  satisfy 


Aw  =  Aw , 

(22) 

v*A  =  Av*, 

(23) 

(w,v)  =  v*w  =  1  . 

(24) 

Now  consider  a  perturbed  matrix  A  -f-  dA,  where  dA  is  a  matrix  of 
small  perturbations  daij.  The  eigenvalues  and  eigenvectors  of  the 
new  matrix  satisfy 

(A  +  dA)(w  -h  dw)  =  (A  -h  dA)(w  -h  dw) .  (25) 

Expanding  the  products  and  eliminating  second-order  terms  yields 
Aw  -h  A(dw)  -f-  (dA)w  =  Aw  -h  A(dw)  -h  (dA)w ,  (26) 

which  simplifies  to 

A(dw)  -f  (dA)w  =  A(dw)  +  (dA)w .  (27) 

Multiplying  both  sides  by  v*  yields 

v*A(dw)  -h  v*(dA)w  =  Av*(dw)  -1-  (dA)v*w .  (28) 

The  first  term  on  the  left-hand  side  is  the  same  as  the  first  term 
on  the  right-hand  side  (because  of  eq.  23),  and  the  last  term  on 
the  right-hand  side  simplifies  to  dA  (because  of  eq.  24),  leaving 

v*dAw  =  dA.  (29) 

If  dA  contains  only  a  single  nonzero  element  daij,  a  change  in 
only  the  zjth  element  of  A,  we  obtain  the  fundamental  sensitivity 


34 


Caswell 


equation: 


dX 

da, 


—  ViW 


3' 


(30) 


'13 


(The  bar  over  Vi  has  been  ignored  in  most  presentations  of  this 
formula,  including  mine.  It  is  irrelevant  to  the  case  of  the  dominant 
eigenvalue  of  a  population-projection  matrix,  which  always  has  real 
eigenvectors,  but  it  must  be  included  for  calculations  involving 
other  eigenvalues.) 

Equation  (30)  says  that  the  sensitivity  of  A  to  changes  in  is 
proportional  to  the  product  of  the  reproductive  value  of  stage  i 
and  the  representation  of  stage  j  in  the  stable  stage  distribution. 

The  sensitivity  of  A  to  changes  in  other  parameters  can  be  cal¬ 
culated  using  the  chain  rule:  for  some  parameter  x, 


dx 


=E 


d\  da, 


13 


daij  dx 


(31) 


The  sensitivity  of  A  gives  the  effect  of  a  small  additive  change  in 
one  of  the  vital  rates.  The  effect  of  a  small  proportional  change  in 
a  vital  rate  is  given  by  the  elasticity  of  A: 


^  aij  dX 
X  daij 


(32) 


In  addition  to  giving  the  proportional  change  in  A  resulting  from 
a  proportional  change  in  the  aij,  the  elasticities  also  measure  the 
contribution  of  the  aij  to  overall  population  growth  rate.  To  be 
precise,  Y.ij  ^ij  —  1  (for  a  simple  proof,  see  Mesterton-Gibbons 
1993),  and  eij  can  be  interpreted  as  the  proportion  of  A  contributed 
by  aij. 

Elasticities  to  other  parameters  can  also  be  calculated: 


e(x) 


X  dX 
X  dx 

X 

X 


X  5A  daij 


daij 


dx 


(33) 


The  elasticities  of  A  with  respect  to  other  parameters  do  not  in 
general  sum  to  one,  and  they  cannot  be  interpreted  as  contributions 
to  population  growth  rate. 
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Sensitivity  or  Elasticity? 

Some  authors  sociu  to  believe  that  sensitivities  and  elasticities  are 
alternatives,  that  one  is  superior  to  the  other,  or  that  one  or  the 
other  is  biased  in  some  way.  This  is  not  so;  they  provide  accurate 
answers  to  different  questions.  The  difference  between  them  is  com¬ 
parable  to  the  difference  between  plotting  the  same  set  of  numbers 
on  arithmetic  (sensitivity)  or  logarithmic  (elasticity)  axes.  Neither 
kind  of  graph  is  wrong,  but  one  or  the  other  may  be  better  at  re¬ 
vealing  interesting  patterns  in  the  numbers.  For  more  discussion, 
see  Chapter  7,  by  Horvitz  et  ai. 


Life-Table-Response  Experiments  and  Comparative  Demography 

Life-table-response  experiments  (LTRE’s)  are  manipulative  exper¬ 
iments  or  comparative  observations  in  which  the  dependent  vari¬ 
able  is  a  complete  set  of  vital  rates  (loosely  speaking,  a  life  ta¬ 
ble;  Caswell  19896).  The  different  environmental  conditions  (the 
‘‘treatments”)  cause  changes  in  the  vital  rates,  which  in  turn  af¬ 
fect  population  dynamics.  LTRE’s  are  often  summarized  by  us¬ 
ing  the  rate  of  increase,  A,  as  a  demographic  statistic  to  integrate 
the  treatment  effects  on  survival  and  reproduction  throughout  the 
life  cycle. 

Knowing  that  a  treatment  produces  a  particular  value  of  A  leaves 
unresolved  the  question  of  how  the  manifold  changes  in  the  vital 
rates  contribute  to  the  effect  on  growth  rate.  After  all,  some  vital 
rates  can  be  changed  a  great  deal  without  affecting  A  (e.g.,  the 
survival  of  a  post-reproductive  age  class),  whereas  small  changes 
in  other  vital  rates  produce  large  changes  in  A.  In  addition,  most 
environmental  factors  have  differential  effects  on  the  different  vital 
rates.  A  given  treatment  may  affect  survival,  growth,  and  fertility 
differently,  with  different  effects  on  those  rates  in  different  stages. 

Treatment  effects  on  A  can  be  decomposed  into  contributions 
from  the  effects  on  each  of  the  vital  rates  (Caswell  19896,  1996u). 
This  decomposition  makes  it  possible  to  pinpoint  where  in  the 
life  history  the  treatment  has  its  greatest  impact.  The  decom¬ 
position  uses  a  first-order  linear  approximation  to  the  effect  on 
A.  I  outline  the  simplest  case  here:  a  set  of  M  treatments  Tm^ 
m  =  1, ...  5  M ,  each  of  which  produces  its  own  matrix  Am  3nd 
population  growth  rate  (I  use  parenthetical  superscripts  to 

denote  treatments  when  subscripts  distinguish  matrix  elements.) 
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Choose  some  condition  as  a  reference  treatment:  this  might  be  the 
mean  of  all  the  treatments,  a  control  treatment,  or  any  other  con¬ 
dition  of  particular  interest.  The  reference  matrix  is  denoted  Ar, 
and  treatment  effects  on  A  are  measured  relative  to 
To  first  order,  we  can  write 


(34) 


(Am  +Ar)/2 


for  m  =  1, . . . ,  M.  Each  term  in  the  summation  gives  the  contri¬ 
bution  of  the  effect  of  treatment  m  on  one  of  the  vital  rates  to  its 
overall  effect  on  A.  That  contribution  is  the  product  of  the  vital- 
rate  effect  and  the  sensitivity  of  A  to  changes  in  that  vital  rate.  If 
either  of  these  terms  is  small — if  the  treatment  doesn’t  effect  aij 
or  if  A  is  insensitive  to  aij — then  the  contribution  of  effects  on 
to  effects  on  A  is  small  The  converse  is  also  true. 

The  sensitivities  in  (34)  must  be  calculated  from  some  particular 
matrix;  here  they  are  calculated  from  a  matrix  “halfway  between” 
the  two  matrices  (A^  and  Ar)  being  compared.  There  is  some 
theoretical  justification  for  this  (Caswell  19896),  and  it  works  well 
in  practice. 

The  reason  for  using  A  as  a  statistic  to  summarize  the  results  of 
an  LTRE  is  that  it  integrates  the  diverse  and  stage-specific  effects 
of  the  treatments.  The  decomposition  analysis  complements  this 
use:  it  pinpoints  the  source,  within  the  life  cycle,  of  the  effects  on 
A.  Experience  with  this  kind  of  analysis  shows  that  it  is  not  safe  to 
assume  that  the  biggest  changes  in  vital  rates  are  responsible  for 
the  effects  of  a  treatment  on  A.  Without  some  analysis  like  (34), 
half  of  the  information  contained  in  an  LTRE  is  wasted. 

Equation  (34)  describes  a  simple,  one-way,  fixed-effect  experi¬ 
mental  design.  The  approach  has  been  extended  to  factorial  de¬ 
signs,  random  designs,  and  regression  designs  (Caswell,  in  press). 
It  can  also  be  applied  to  statistics  other  than  A  (as  long  as  a  per¬ 
turbation  theory  is  available  for  the  statistic)  and  to  parameters 
other  than  matrix  elements  (Caswell  1989c,  1996a,  in  press). 


Prospective  and  Retrospective  Analyses 

The  preceding  subsections  outline  two  ways  of  using  perturbation 
analysis.  Sensitivity  and  elasticity  calculations  are  prospective  anal¬ 
yses;  they  predict  the  results  of  perturbations  of  the  vital  rates 
before  they  happen.  Indeed,  they  even  show  the  results  of  pertur¬ 
bations  that  are  biologically  impossible.  They  tell  nothing  about 
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which  vital  rates  are  actually  responsible  for  an  observed  change 
in  A.  The  LTRE  decomposition  analysis  answers  this  kind  of  retro¬ 
spective  question.  It  does  so  by  combining  sensitivity  analysis  with 
information  on  the  actual  variance  in  the  aij. 

Don’t  confuse  these  two  kinds  of  analysis,  especially  in  ambigu¬ 
ous  questions  such  as,  ‘‘which  of  the  vital  rates  is  most  important 
to  population  growth?”  One  way  to  answer  this  question  is  to  find 
the  rate  with  the  biggest  sensitivity  (or  elasticity);  that  rate  is  most 
important  in  the  sense  that  if  you  were  to  change  all  the  rates  by 
the  same  amount  (or  same  proportion),  it  would  have  the  biggest 
impact-  Another  answer  is  based  on  an  LTRE  analysis;  the  most 
important  vital  rate  is  the  one  with  the  variation  that  makes  the 
biggest  contribution  to  the  variability  in  A.  The  two  answers  are 
usually  different-  Both  are  valid,  but  they  answer  different  ques¬ 
tions,  the  first  prospective,  the  second  retrospective.  Chapter  7,  by 
Horvitz  et  ah,  explores  these  issues  further. 


4  Density-Dependent  Matrix  Models 

The  models  analyzed  so  far  have  been  linear  and  time-invariant; 
the  vital  rates  axe  independent  of  population  density  or  tempo¬ 
ral  changes.  Time- varying  models  are  discussed  in  Chapter  3;  here 
I  consider  the  inclusion  of  density  dependence,  which  makes  the 
model  nonlinear.  Although  the  mathematical  tools  are  different 
from  those  used  in  the  linear  case,  the  focus  is  still  on  asymp¬ 
totic  behavior  and  perturbation.  Unlike  linear  models,  density- 
dependent  models  do  not  grow  exponentially-  Instead,  solutions 
tend  to  converge  to  limited  subsets  of  the  state  space,  called  at¬ 
tractors.  The  attractors  may  he  fixed  points  (also  called  equilibria)^ 
cycles,  or  more-complicated  structures. 

Two  kinds  of  perturbation  analysis  are  important.  One  asks  the 
effect  of  perturbing  initial  conditions.  This  is  of  no  interest  in  the 
linear  case,  since  the  ergodic  theorem  guarantees  convergence  to 
the  stable  population  structure  and  asymptotic  growth  rate  from 
any  nonzero  initial  condition.  In  a  density-dependent  model,  how¬ 
ever,  small  perturbations  of  initial  conditions  can  lead  to  very  dif¬ 
ferent  dynamics,  depending  on  the  stability  of  a  fixed  point.  The 
second  kind  of  perturbation  considers  changes  in  the  parameters  of 
the  model.  Often,  small  changes  in  parameters  leave  the  qualitative 
asymptotic  behavior  unchanged.  But  sometimes  small  parameter 
changes  have  big  effects  on  asymptotic  dynamics:  stable  trajecto¬ 
ries  become  unstable,  fixed  points  are  created  or  destroyed,  attrac- 
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tors  change  from  fixed  points  to  cycles,  cycles  give  way  to  chaos, 
etc.  These  qualitative  changes  are  called  bifurcations,  and  finding 
them  is  a  major  pastime  of  people  who  study  density-dependent 
models. 

There  are  many  ways  of  incorporating  density  dependence  in  ma¬ 
trix  models  (Caswell  1989a;  see  also  Cushing  1988;  Getz  &  Haight 
1989;  Silva  &  Hallam  1992;  Logofet  1993;  Dennis  et  al.  1995).  “Den¬ 
sity”  can  be  defined  as  the  abundance  of  a  stage  or  set  of  stages, 
as  a  weighted  combination  of  the  abundances  of  a  set  of  stages,  or 
simply  as  the  total  number  ^  .  Density  can  affect  the  vital 

rates  at  many  points  in  the  life  cycle  or  at  only  one.  Each  of  the 
vital  rates  may  be  a  different  function  of  the  density,  or  the  entire 
set  of  vital  rates  may  be  affected  by  density  in  the  same  way.  In  this 
section,  I  introduce  three  models  that  incorporate  simple  density 
effects  on  reproduction,  on  growth,  and  on  survival,  using  them  to 
demonstrate  some  dynamic  consequences  of  density  dependence. 


Basic  Formulations 

Consider  a  two-stage  model,  with  a  matrix 


Pi  F2 

Gi  P2 


and  suppose  that  the  matrix  entries  are  given  by 


Pi  =  C7i(l*-7i), 

Gi  =  cri7i, 

Po  =  G  2  . 


(35) 


where  <j\  and  <72  are  survival  probabilities,  and  71  is  the  probabiliry 
of  growth  from  stage  1  to  stage  2,  approximated  by  1/ri.  where  Ti 
is  the  mean  duration  of  stage  1. 

This  simple  model  contains  the  rates  of  reproduction,  growth, 
and  survival,  each  of  which  can  be  made  density-dependent. 


Density- dependent  reproduction.  Let  fertility  at  zero  density  be 
given  by  /o,  and  suppose  that  fertility  declines  exponentially  with 
increases  in  density.  Then, 


F2{N)  =  foexp{-bN), 


(36) 
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where  6  is  a  constant  that  measures  the  strength  of  density  depen¬ 
dence.  The  resulting  matrix  is 


Pi  foe-^^  ■ 

Gi  P2  / 


(37) 


Density-dependent  survival  Let  0*1^0  sind  ^2,0  denote  the  survival 
probabilities  of  stages  1  and  2  at  zero  density.  Assume  that  both 
juvenile  and  adult  survival  probabilities  are  affected  in  the  same 
way  by  density,  so  that 

ai{N)  =  ai,oexp(-6iV) ,  (38) 

0-2  (A)  =  (72,0  exp(-6iV) .  (39) 

The  resulting  matrix  is 


0-1,0(1-71)6  F2 

0-1,0716"''^  cr2,oe~^^ 

■  F2 

Gi,oe-’>^^  P2,oe-^^  ' 


(40) 


Density- dependent  growth.  Suppose  that  the  mean  duration  of 
the  juvenile  stage  increases  with  density: 

r[N)  =  roexp(feiV); 

then  the  growth  probability  is  given  by 


l{N)  = 


1 

t{N) 

7o  exp{-hN). 


(41) 


The  resulting  matrix  is 


A(g) 


<7i(l  -7oe  *’^) 

F2 

P2 

0-1(1  -  7oe~'’^ 

F2 

P2 

(42) 


The  assumption  of  exponential  density  dependence  (“overcompen¬ 
sation”  in  the  language  of  fisheries  biology)  has  important  dynamic 
consequences.  A  strictly  compensatory  form  (e.g.,  1/(1  +  hN)).,  or 
a  depensatory  form  with  some  range  of  positive  density  depen¬ 
dence,  may  produce  different  dynamics  (Caswell  1989a;  Silva  & 
Haliam  1992), 
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Equilibria  and  Stability 

The  equilibria  of  a  density-dependent  model  are  population  vectors 
fi  that  satisfy 

n  =  A^n.  (43) 

Depending  on  the  form  of  the  density  dependence,  there  may  be 
more  than  one  equilibrium.  In  the  absence  of  immigration,  n  = 
0  is  always  an  equilibrium.  With  a  little  luck,  equation  (43)  can 
be  solved  analytically  for  fi.  More  often  than  not,  however,  the 
equilibria  must  be  found  numerically. 


Density-dependent  reproduction.  The  equilibrium  n  is  defined  by 


fii  =  Pini +/o€‘‘'^’n2 : 
n2  —  GiUi  +  P2n2  . 

The  second  of  these  equations  can  be  solved  for  ni: 


ni  = 


I-P2. 


"712; 


substituting  this  in  the  first  equation  yields 

Combining  these  two  results  gives  an  expression  for  h2, 

GiN 

”2"  (I  +  G1-P2)’ 
and  ni  can  be  found  as  fii  =  iV  —  n2. 


(44) 

(45) 

(46) 

(47) 

(48) 


(49) 


Density-dependent  survival  The  equations  defining  the  equilib¬ 
rium  are 

fii  =  Pi,oe“^^ni  -h  F2n2  , 

712  =  +  P2,0^”^'^fi2  • 

The  second  of  these  equations  can  be  solved  for  722: 

..  Gi^Qe~~^^  ^  ^  ^ 

(l-P2,o)e-^^  ' 

When  substituted  into  the  first  equation  this  eventually  yields 

Pl,0  ““  Pl,oP2,0^'”^^  +  P2Gi,0  +  P2,0  =  •  (51) 
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Multiplying  both  sides  by  gives  a  quadratic  equation  in  the 
new  variable  y  =  ,  Solving  this  equation  for  y  and  substituting 

leads  to 

N-logyb, 

=  +  ,  (52) 

n2  =  N  —  hi. 


Density- dependent  growth.  The  equations  to  be  solved  for  the 
equilibrium  are 


hi  -ai(l-  7oe  j  hi  +  F2^2  , 
h2  =  cri7oe“^^ni  -h  P2^2  • 


(53) 


As  before,  the  second  equation  can  be  solved  for  n2, 


n2 


<Ti7oe-^^ 

I-P2 


ni ; 


(54) 


substituting  this  relation  into  the  first  equation  eventually  leads  to 


V^i7o(F2  +  P2  1)/ 


(55) 


None  of  these  analytical  solutions  for  n  is  particularly  informa¬ 
tive  at  first  glance.  This  is  typical;  only  in  exceptional  circum¬ 
stances  are  the  formulas  for  the  equilibria  in  a  matrix  population 
model  simple  enough  to  appear  informative. 

An  equilibrium  is  said  to  be  locally  stable  if  small  perturbations 
remain  close  to  the  equilibrium,  and  locally  asymptotically  stable 
if  small  perturbations  eventually  return  to  the  equilibrium.  The 
adjective  “local”  refers  to  the  smallness  of  the  perturbations;  it  may 
well  happen  that  small  perturbations  return  to  an  equilibrium,  but 
large  ones  are  attracted  to  another  equilibrium,  or  to  some  other 
kind  of  attractor  (a  cycle,  for  example,  or  a  strange  attractor). 

The  local  stability  of  a  fixed  point  is  determined  by  approximat¬ 
ing  the  nonlinear  density-dependent  model  by  a  linear  model  that 
is  accurate  for  small  perturbations.  Begin  by  defining  a  vector  x  of 
deviations  from  the  equilibrium  n: 


x(t)  =  n{t)  -  n . 


(56) 
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The  dynamics  of  x,  near  the  equilibrium,  are  given  by 

x(t  +  l)  =Bx(t).  (57) 

where  the  constant  matrix  B,  called  the  Jacobian  matrix,  is  the 
linear  approximation  to  the  full  nonlinear  system  near  the  equi¬ 
librium. 

To  illustrate,  let  us  write  a  two-dimensional  matrix  model  as 


ni(t-l-l)  =  uiini(t) -(-ai2n2(t)  = /i(n) ,  (58) 

^2(^  +  1)  =  o-2ini{t)  +  a22n2{t)  =  f2{n) ,  (59) 


with  the  understanding  that  aij  =  ay(n). 
The  Jacobian  matrix  B  is 


B  = 


dni 

dh 

dni 


Ml 

dn2 

ih, 

dn2 


Differentiating  /i  gives 
dni 

dn2 


011  4-  ni 


da 


dl 


dni 


+  n2 


da 


12 


dan 

—  ai2  -I-  ni  — - h  n2 


dni 

dai2 


dn2  '  dn2 

The  derivatives  of  /2  have  the  same  form.  Thus,  B  is  given  by 

da 


(60) 

(61) 

(62) 


B  —  Aft  -H 


=  A"  +  [  ^“  •  (63) 

All  of  the  derivatives  are  evaluated  at  the  equilibrium. 

Equation  (63)  makes  the  calculation  of  the  Jacobian  straightfor¬ 
ward  for  matrix  population  models  (it  is  due  to  Beddington  1974). 
Given  the  equilibrium  vector  fi,  calculate  the  derivatives  of  A  with 
respect  to  each  of  the  ui  and  multiply  these  matrices  by  n;  put 
the  resulting  vectors  as  columns  in  a  matrix,  and  add  this  to  A 
evaluated  at  the  equilibrium. 

The  equilibrium  n  is  asymptotically  stable  if  the  eigenvalues  of 
B  are  all  less  than  one  in  magnitude  (or  are  '^within  the  unit  cir¬ 
cle,”  referring  to  the  circle  with  radius  one  in  the  complex  plane). 
It  is  imstable  if  any  of  the  eigenvalues  are  outside  the  unit  circle. 
An  eigenvalue  falling  exactly  on  the  unit  circle  (i.e.,  with  magni¬ 
tude  exactly  equal  to  one)  signals  a  transition  from  stability  to 
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instability,  or  vice  versa,  and  is  considered  in  the  next  section. 

Remember  that  (63)  describes  the  dynamics  of  perturbations, 
not  the  dynamics  of  numbers  of  individuals.  Thus,  the  entries  of 
X  can  be  negative,  and  the  matrix  B  often  contains  negative  en¬ 
tries.  The  Perron-Ftobenius  theorem  is  of  little  use  in  evaluating 
the  eigenvalues  of  B,  except  in  special  cases  (DeAngelis  et  al.  1980; 
Caswell  1989a,  Example  9.3),  and  the  eigenvalue  of  largest  magni¬ 
tude  may  be  negative  or  complex. 

The  Jacobian  matrices  for  our  three  density-dependent  models 
can  be  written  down  easily.  Note  that  all  density  effects  depend  on 
AT  =  ni  +  n2,  so  that 

^o,ij  Oclij  ddij  f64') 

dni  dn2  dN 

Thus,  the  matrix  B  reduces  to 

+  "■].  (65) 

oN  I  Tl2  712  _ 

The  resulting  Jacobian  matrices  for  the  three  example  models  are 
as  follows:  for  density-dependent  reproduction, 

for  density-dependent  survival, 


Pi,oe-^^  F2 

Gx,oe-^^  P2,oe-^^  _ 

f^-bN  Pl,oni  Pl,0^1 

+  P2,0^2  +  P2,0^2 


and  for  density-dependent  growth, 

[  cri(l  -706’'’^)  F2 
Gi,oe-^^  P2 


ni  ni 

—Til  —ni 


In  each  case,  N  and  n  are  evaluated  using  the  appropriate  equi¬ 
librium  formulas.  Note  that  I  have  not  attempted  to  insert  the 
formulas  for  the  equilibria.  In  general,  it  is  a  lucky  circumstance  in 
which  such  expressions  simplify  usefully.  More  often,  the  best  that 
can  be  done  is  to  write  the  Jacobian  in  terms  of  the  equilibria  and 
to  carry  out  further  analyses  numerically. 
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In  studying  a  linear  model,  the  response  of  the  eigenvalues  to 
changes  in  parameters  is  usually  more  interesting  than  their  values 
for  one  specific  parameter  set.  Similarly,  in  studying  a  nonlinear 
model  it  is  usually  more  interesting  to  see  how  the  stability  proper¬ 
ties  change  as  parameters  are  varied  than  to  characterize  stability 
for  one  parameter  set.  A  qualitative  change  in  stability  is  called  a 
bifurcation. 


Bifurcations 

Roger  Tory  Peterson  published  the  first  edition  of  A  Field  Guide  to 
the  Birds  in  1934.  At  the  time,  the  standard  guide  for  the  would-be 
birdwatcher  was  Chapman’s  (1932)  Handbook  of  Birds  of  Eastern 
North  America,  first  published  in  1895.  Chapman’s  book  contained 
detailed  dichotomous  keys  to  the  bird  species,  most  of  which  would 
be  useless  without  having  the  bird  in  one  hand  (and,  most  likely, 
a  smoking  shotgun  in  the  other).  Peterson’s  innovation  was  to  sac¬ 
rifice  rigor  but  to  focus  on  characters  that  would  distinguish  most 
of  the  species,  most  of  the  time,  in  the  field. 

My  goal  in  this  section  is  to  provide  a  kind  of  field  guide  to 
bifurcations  of  equilibria  of  nonlinear  matrix  models.  I  describe 
the  basic  ideas  of  bifurcation  theory,  without  assuming  that  my 
reader  is  equipped,  with  the  powerful  tools  of  rigorous  bifurcation 
theory,  but  assuming  the  ability  to  conduct  numerical  simulations. 
The  question  is,  what  to  look  for  that  will  help  to  identify  the 
bifurcations  that  occur. 

Be  forewarned  that,  just  as  there  are  situations  where  Peterson 
had  to  admit  that  a  reliance  on  field  marks  was  misleading  or 
impossible,  there  are  no  doubt  situations  in  which  my  suggestions 
here  will  be  incorrect.  (In  a  notorious  case,  Chapman  distinguished 
the  Alder  and  Least  flycatchers  of  the  genus  Empidonax  on  i  the 
basis  of  whether  the  wing  is  longer  or  shorter  than  2.60  inches 
and  whether  the  lower  mandible  is  flesh-colored  or  strongly  tinged 
with  brownish.  Peterson,  by  contrast,  said  that  “it  is  quite  risky  to 
attempt  to  tell  them  apart  by  mere  variations  in  color”  and  relied 
on  calls  during  the  breeding  season.)  In  such  cases,  you  can  make 
use  of  the  many  detailed  treatments  of  bifurcation  theory  (two 
excellent  recent  examples  are  Wiggins  1990  and  Hale  &  Kogak 
1991;  see  also  the  review  paper  Whitley  1983)  that  provide  the 
mathamatical  analogue  of  the  birdwatcher’s  shotgun.  See  Chapter 
6,  by  Cushing,  for  an  introduction  to  the  application  of  this  theory 
to  matrix  models. 
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Types  of  bifurcations.  An  equilibrium  point  can  bifurcate  when  it 
loses  stability  in  response  to  a  change  in  a  parameter,  that  is,  when 
the  dominant  eigenvalue  A  of  the  Jacobian  matrix  B  in  equation 
(63)  is  equal  to  one.  This  can  happen  in  three  ways: 

A  =  1,  a  “+1  bifurcation,” 

A  =  -1,  a  bifurcation,”  (67) 

A  =  a  ±  a  ‘‘complex  conjugate  pair”  bifurcation. 

When  an  eigenvalue  (or  a  complex  conjugate  pair)  crosses  the 
unit  circle,  the  others  must  still  be  inside,  or  there  would  be  no 
change  in  stability,  since  the  equilibrium  would  already  be  unsta¬ 
ble.  (I  am  ignoring  the  rare  cases  in  which  more  than  one  eigen¬ 
value  crosses  the  unit  circle  simultaneously.)  Perhaps  surprisingly, 
the  bifurcation  does  not  depend  on  what  those  other  eigenvalues 
are  doing,  nor  on  how’  many  of  them  there  are.  The  mathematical 
expression  of  this  fact  is  the  Center  Manifold  Theorem  (e.g.,  Wig¬ 
gins  1990).  Roughly,  it  says  that  the  state  space  can  be  divided 
into  two  parts,  one  associated  with  the  eigenvalue  leaving  the  unit 
circle  and  the  second  associated  with  all  the  other  eigenvalues.  The 
dynamics  on  the  second  part  remain  stable,  while  the  bifurcation 
can  be  studied  on  the  first  part.  Thus  all  the  basic  behaviors  of 
-fl  and  —1  bifurcations  can  be  studied  in  one-dimensional  mod¬ 
els,  while  the  Hopf  bifurcation  can  be  studied  in  two-dimensional 
models. 

What  different  bifurcations  are  possible  in  response  to  a  param¬ 
eter  change?  First,  H-l  bifurcations  are  of  three  types. 

1.  At  a  transcritical  bifurcation,  two  fixed  points,  one  stable  and 
the  other  unstable,  collide  and  exchange  stability.  This  occurs  in 
all  density-dependent  matrix  models  as  a  bifurcation  of  the  zero 
equilibrium.  On  one  side  of  the  bifurcation  the  zero  equilibrium 
is  stable.  An  unstable  equilibrium  exists,  but  it  is  negative  and 
usually  ignored  in  our  thinking  about  populations.  However,  it 
does  exist.  On  the  other  side  of  the  bifurcation,  the  zero  equi¬ 
librium  is  unstable,  and  there  is  a  stable  positive  equilibrium. 
See  Chapter  6,  by  Cushing,  for  a  detailed  description  of  this 
bifurcation. 

2.  On  one  side  of  the  saddle-node  bifurcation,  there  is  no  fixed 
point  at  all.  At  the  bifurcation  point,  two  fixed  points  appear, 
one  stable  and  the  other  unstable.  What  you  see  in  a  saddle- 
node  bifurcation  depends  on  where  you  look.  If  you  follow  the 
stable  fixed  point,  it  suddenly  disappears  (as  it  collides  with  the 
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unstable  fixed  point  of  which  you  were  blissfully  unaware),  and 
trajectories  move  off  to  some  other  attractor  in  the  state  space. 

Saddle-node  bifurcations  occur  in  population  models  when  there 
is  an  Allee  effect  (also  known  as  depensatory  density  depen¬ 
dence)  present,  that  is,  when,  at  low  densities,  survival  or  fertil¬ 
ity  depends  positively  on  density.  Typically,  such  a  population 
has  three  fixed  points:  zero,  an  intermediate  critical  density,  and 
a  high  density.  Zero  and  the  high  density  are  both  stable  fixed 
points,  while  the  intermediate  critical  density  is  unstable.  Popu¬ 
lations  below  the  critical  density  decay  to  extinction,  and  those 
above  it  increase  to  the  high-density  fixed  point.  As  some  param¬ 
eter  is  varied,  the  high-density  fixed  point  and  the  intermediate- 
density  fixed  point  collide  and  annihilate  each  other;  the  popu¬ 
lation  at  that  point  crashes  to  extinction.  Cushing  (Chapter  6) 
shows  an  example. 

3.  On  one  side  of  a  pitchfork  bifurcation,  there  is  one  stable  fixed 
point.  On  the  other  side,  one  unstable  fixed  point  is  surroimded 
by  two  new,  stable  fixed  points.  At  the  bifurcation,  the  two  sta¬ 
ble  fixed  points  appear  close  to  the  now-unstable  fixed  point;  a 
trajectory  may  approach  either  of  the  two  stable  fixed  points,  de¬ 
pending  on  initial  conditions.  The  canonical  examples  of  pitch- 
fork  bifurcations  occur  in  one-dimensional  cubic  maps.  These 
maps  appear  in  some  genetic  models  (May  1979),  but  pitchfork 
bifurcations  rarely  appear  in  density-dependent  matrix  models. 

Second  are  -1  bifurcations,  or  flip  bifurcations.  At  this  type 
of  bifurcation  point,  a  stable  fixed  point  becomes  imstable,  and 
a  stable  2-cycle  appears.  Just  beyond  the  bifurcation  point,  the 
amplitude  of  the  2-cycle  is  small,  and  it  surrounds  the  unstable 
fixed  point. 

The  third  type  of  bifurcation  develops  from  complex-conjugate 
pairs  of  eigenvalues  and  is  also  called  a  Hopf  or  Naimark-Sacker 
bifurcation.  As  a  complex  conjugate  pair  of  eigenvalues  leave  the 
imit  circle,  the  fixed  point  loses  its  stability  and  an  invariant  circle 
forms  around  the  now-imstable  fixed  point.  An  invariant  circle  is  a 
continuous  closed  curve  (not  necessarily  an  exact  circle).  Any  point 
on  this  curve  is  mapped  to  another  point  on  the  curve,  hence  the 
term  ‘invariant.” 

Trajectories  on  the  invariant  circle  may  be  periodic,  or  they  may 
be  quasiperiodic,  rotating  around  the  circle  without  ever  repeat¬ 
ing  themselves,  depending  on  the  location  of  the  eigenvalues  when 
they  cross  the  unit  circle  at  the  bifurcation.  In  general,  trajecto- 


MATRIX  METHODS  FOR  POPULATION  ANALYSIS 


47 


lies  immediately  following  a  Hopf  bifurcation  are  quasiperiodic, 
containing  two  periods  that  are  not  integer  multiples  of  each  other 
(see  example  below). 

There  are  two  other  possibilities.  One  is  called  strong  resonance 
and  occurs  when  the  bifurcation  point  happens  at 

A  =  where  k  =  l,  V2,  Vs,  V4.  (68) 

The  second  is  called  weak  resonance  and  occurs  when 

A  =  (69) 

where  p/q  is  a  ratio  of  small  integers  but  does  not  equal  any  of  the 
first  four  roots  of  unity.  In  this  case,  the  trajectory  on  the  invariant 
circle  has  period  q. 

Strong  resonance  is  complicated.  Mathematicians  say  things  like 
“the  strong  resonances  . . .  exhibit  rather  different  behavior,  the 
details  of  which  are  not  yet  fully  understood”  (Whitley  1983),  or 
“the  dynamics  of  such  maps  . . .  can  be  exceedingly  complicated 
and  the  answer  is  not  yet  completely  known”  (Hale  &  Ko^ak  1991). 
I  cannot  improve  on  that,  but  it  appears  that  bifurcations  with 
strong  resonance  often  have  trajectories  with  period  1/k  lurking  in 
their  dynamics,  often  mixed  up  with  quasiperiodic  and,  eventually, 
chaotic  dynamics  (for  an  example,  see  Guckenheimer  et  al.  1977). 

Examples:  bifurcations  of  the  zero  equilibrium.  The  three  simple 
density-dependent  models  in  this  section  provide  examples  of  these 
bifurcations.  Zero  is  an  equilibrium  of  all  three  models;  to  examine 
its  bifurcations,  evaluate  the  Jacobian  matrix  at  zero  and  track  the 
magnitude  of  its  eigenvalues  as  parameters  in  the  matrix  change.  In 
the  examples  below,  I  show  bifurcations  resulting  from  changes  in 
reproductive  output  F2  (or  /o  in  the  model  with  density-dependent 
fertility).  Note,  however,  that  Cushing  (1988;  Chapter  6)  has  shown 
in  general  that  the  bifurcation  parameter  can  be  taken  to  be  the 
net  reproductive  rate  (the  expected  number  of  offspring  produced 
by  an  individual  over  its  lifetime);  he  has  shown  how  to  calculate 
this  quantity  for  arbitrary  stage-structured  matrices. 

The  Jacobian  matrix  for  the  zero  equilibrium  is  particularly  sim¬ 
ple,  because  the  second  term  on  the  right-hand  side  of  equation 
(63)  is  zero.  Thus,  the  Jacobian  is  simply  the  projection  matrix 
evaluated  at  zero,  and  it  is  the  same  in  all  three  examples. 

Figure  3  shows  the  magnitude  of  the  dominant  eigenvalue  of 
B  as  a  function  of  F2]  it  exceeds  one  when  F2  ^  0.8.  The  lower 
panel  shows  the  location  in  the  complex  plane  of  the  dominant 
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Figure  3.  The  bifurcation  of  the  zero  equilibrium  for  the  three  two- 
stage  examples.  Upper  panel,  The  magnitude  of  the  dominant  eigenvalue 
of  the  Jacobian  matrix  B  as  a  function  of  the  reproductive  output  F2. 
Lower  panel.  The  location,  relative  to  the  unit  circle  in  the  complex 
plane,  of  the  eigenvalue  just  before  (open  circle)  and  after  (closed  circle) 
the  bifurcation.  Parameter  values:  ax  =  0.6,  ^2  =  0.8,  71,0  =  0.2,  6  =  1. 

eigenvalue  before  and  after  the  bifurcation;  this  is  a  +1  bifurcation 
since  the  eigenvalue  leaves  the  unit  circle  at  -f  1.  Figure  4  shows  the 
equilibrium  (for  total  population  size  iV  =  m  +  712)  as  a  function 
of  F2.  The  pattern  is  typical  of  a  transcriticai  bifurcation.  Note 
that  all  three  models  bifurcate  from  zero  at  the  same  value  of  ^2, 
because  all  three  matrices  are  identical  evaluated  at  n  =  0,  even 
though  they  differ  in  their  subsequent  behavior. 

Bifurcations  of  the  positive  equilibrium.  What  happens  to  the  sta¬ 
ble  positive  equilibrium  that  appears  after  the  zero  equilibrium 
becomes  imstable?  Each  of  the  three  models  shows  a  different  bi¬ 
furcation  pattern.  I  show  examples  for  each  model,  using  F2  (or  /o 
in  the  model  with  density-dependent  fertility)  as  the  bifurcation 
parameter.  There  is  no  reason  to  expect  that  this  exhausts  the 
possibilities;  other  types  of  bifurcations  may  well  be  produced  by 
varying  other  parameters. 

For  density-dependent  reproduction,  the  bifurcation  occurs  when 
/o  95.  The  eigenvalue  leaves  the  unit  circle  at  -1  (Fig.  5), 
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Figure  4.  Bifurcation  plots  for  the  three  two- stage  examples,  showing  the 
equilibrium  value  of  N  =  ni  +  n2  as  a  function  of  reproductive  output 
(F2;  in  the  model  with  density- dependent  fertility,  fo)-  Other  parameter 
values  as  in  Figure  3, 
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Figure  5.  The  bifurcation  of  the  positive  equilibrium  for  the  density- 
dependent  fertility  model  as  a  function  of  reproductive  output  fo  -  Inter¬ 
pretation  and  parameter  values  as  in  Figure  3. 
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Figure  6.  Bifurcation  plots  for  the  three  two-stage  examples,  showing  the 
long-term  behavior  of  total  population  size  (N  =  ni+n2^  as  a  function  of 
reproductive  output  (F2;  for  the  model  with  density- dependent  fertility, 
fo ).  Other  parameters  as  in  Figure  3. 

implying  a  flip  bifurcation,  as  seen  in  Figure  6a.  The  fixed  point  is 
replaced  by  a  2-cycle. 

For  the  model  with  density-dependent  survival,  the  equilibrium 
becomes  unstable  when  F2  ~  50.  The  eigenvalue  leaves  the  unit  cir¬ 
cle  very  near  (but  not  quite  on)  the  imaginary  axis  (Fig.  7).  This 
implies  a  Hopf  bifurcation,  but  one  very  close  to  strong  resonance. 
The  bifurcation  pattern  is  shown  in  Figure  66.  Just  after  the  bifur¬ 
cation  point,  trajectories  are  quasiperiodic  on  an  invariant  circle. 
As  F2  increases  further,  the  invariant  circle  begins  to  look  more 
like  an  invariant  square,  finally  collapsing  to  a  stable  4-cycle  (Fig. 
8).  See  Wikan  and  Mjolhus  (1995)  for  an  exploration  of  this  period- 
4  behavior  in  a  fully  age-classified  model  with  density-dependent 
survival. 

For  the  model  with  density-dependent  growth,  the  bifurcation 
occurs  when  F2  «  1200.  This  time,  the  eigenvalues  leave  the  unit 
circle  as  a  complex  conjugate  pair  (Fig.  9)  that  is  not  one  of  the  first 
four  roots  of  unity.  The  result  is  a  Hopf  bifurcation  to  quasiperiodic 
dynamics  on  an  invariant  circle  (Fig.  10). 

Subcritical  and  supercritical  bifurcations.  The  flip,  pitchfork,  and 
Hopf  bifurcations  come  in  two  varieties,  the  supercritical  (which 
I  have  discussed  so  far)  and  the  subcritical.  For  purposes  of  this 
discussion,  suppose  that  the  bifurcation  happens  when  /  =  /c,  and 
that  the  fixed  point  is  stable  for  f  <  fc  and  unstable  for  /  >  /c.  In 
the  supercritical  bifurcations,  the  imstable  fixed  point  for  f  >  fc 
is  accompanied  by  a  stable  2-cycle,  two  stable  fixed  points,  or  a 
stable  invariant  circle.  In  the  subcritical  bifurcations,  the  2-cycle, 
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Figure  7.  The  bifurcation  of  the  positive  equilibrium  for  the  model  with 
density-dependent  survival  as  a  function  of  reproductive  output  Fa.  In¬ 
terpretation  and  parameter  values  as  in  Figure  3. 


Figure  8.  Phase  portraits  for  the  asymptotic  dynamics  of  the  model  with 
density-dependent  survival  The  inner  invariant  curve  is  for  F2  =  52; 
the  outer  curve,  F2  =  53;  the  Fa  =  55.  Other  parameters  as  in 

Figure  3. 
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Figure  9.  The  bifurcation  of  the  positive  equilibrium  for  the  model  with 
density -dependent  growth  as  a  function  of  reproductive  output  Fb.  Inter¬ 
pretation  and  parameter  values  as  in  Figure  3. 


Figure  10.  Phase  portraits  for  the  asymptotic  dynamics  of  the  model 
with  density -dependent  growth.  The  inner  equilibrium  point  is  for  F2  = 
1150;  the  invariant  curves  are  for  F2  =  1250  and  F2  —  1300.  Other 
parameters  as  in  Figure  3. 
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the  two  fixed  points,  or  the  invariant  circle  exist  for  /  <  /c,  but 
they  are  unstable.  At  the  bifurcation  point  /c,  the  original  fixed 
point  becomes  unstable  and  the  other  structures  disappear. 

What  happens  after  a  subcritical  bifurcation  (i.e.,  for  /  >  /c) 
is  hard  to  predict.  Since  the  fixed  point  becomes  unstable  at  /c, 
trajectories  do  not  remain  close  by,  but  no  nearby  stable  fixed 
points,  cycles,  or  invariant  circles  appear  to  attract  trajectories. 
Therefore,  solutions  leave  the  vicinity  of  the  original  fixed  point, 
and  their  dynamics  depend  on  the  location  of  other  attractors  in 
the  system. 

Subcritical  bifurcations  are  not  as  common  as  supercritical  ones, 
but  they  cannot  be  ignored.  Neubert  and  Kot  (1992)  found  that,  in 
four  discrete  predator-prey  models  in  which  the  prey  alone  exhib¬ 
ited  a  supercritical  flip  bifurcation,  the  system  with  both  predator 
and  prey  exhibited  a  subcritical  flip  bifurcation,  leading  to  ex¬ 
tinction  of  the  predator.  Guckenheimer  et  al.  (1977)  and  Levin 
(1981)  found  subcritical  bifurcations  in  age-structured  models  with 
strongly  resonant  Hopf  bifurcations,  and  Dennis  et  al.  (1995)  found 
a  subcritical  bifurcation  in  a  stage-structured  model  fit  to  labora¬ 
tory  data  on  Tribolium  populations. 

Analytical  methods  exist  for  determining  whether  a  bifurcation 
is  supercritical  or  subcritical,  on  the  basis  of  terms  in  the  Taylor 
series  expansion  of  the  model  near  the  fixed  point.  However,  these 
methods  require  that  the  model  first  be  reduced  to  one  or  two  di¬ 
mensions  using  the  Center  Manifold  Theorem,  which  is  not  always 
easy. 

On  beyond  the  first  bifurcation.  The  positive  fixed  point  of  a  ma¬ 
trix  population  model  usually  first  loses  stability  through  a  super¬ 
critical  flip  or  supercritical  Hopf  bifurcation  and,  then,  is  replaced 
by  a  2-cycle  or  a  quasiperiodic  trajectory  on  an  invariant  circle. 
As  the  bifurcation  parameter  increases  further,  these  attractors  go 
through  a  series  of  bifurcations,  often  leading  to  chaos.  The  flip 
bifurcation  is  typically  followed  by  a  series  of  period-doublings. 
Successively  higher  periods  are  stable  for  smaller  ranges  of  the  bi¬ 
furcation  parameter;  eventually,  trajectories  become  chaotic.  The 
resulting  bifurcation  diagram  is  similar  to  the  familiar  diagram  for 
the  discrete  logistic  equation. 

Quasiperiodic  solutions  produced  by  a  Hopf  bifurcation  usually 
go  through  a  series  of  frequency  locking  events  in  which  trajectories 
become  periodic  for  a  range  of  values  of  the  bifurcation  parameter. 


54 


Caswell 


D-D  Fertility 


D-D  Survival 


D-D  Growth 


Reproductive  Output 


Figure  11.  Bifurcation  plots  beyond  the  first  bifurcation  for  the  three 
two-stage  examples,  showing  the  long-term  behavior  of  total  population 
size  (N  =  ni  n2)  as  a  function  of  reproductive  output  (F2  or,  for 
the  model  with  density- dependent  fertility,  fo)-  Other  parameters  as  in 
Figure  3. 


Frequency  locking  is  followed  by  more  quasiperiodicity.  Eventually, 
the  invariant  circle  may  begin  to  break  up,  or  become  wrinkled  or 
folded,  leading  to  chaotic  dynamics. 

Both  patterns  can  be  found  in  our  examples.  Figure  11  shows 
the  bifurcation  diagrams  for  the  three  models  as  reproductive  out¬ 
put  is  increased  even  further,  and  Figure  12  shows  examples  of 
the  attractors  in  the  chaotic  regime.  Note  that  although  all  three 
models  eventually  become  chaotic,  the  structures  of  the  resulting 
attractors  are  quite  different. 
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Figure  12.  Phase  portraits  for  chaotic  attractors  for  the  three  two- stage 
examples,  plotting  n2  vs.  ni  for  several  thousand  iterations,  after  tran¬ 
sients  have  died  out  For  the  growth  plot,  the  vertical  axis  is  100,  000xti2- 
Values  of  reproductive  output  (F2  or,  for  the  model  with  density- depen¬ 
dent  fertility,  fo )  are  shown;  other  parameters  as  in  Figure  3. 


5  Conclusion 

This  chapter  introduces  some  of  the  most  immediately  useful  meth¬ 
ods  for  analyzing  population  dynamics  with  matrix  models.  In 
combination  with  Chapter  3  on  stochastic  models  (Tuljapurkar), 
and  Chapter  6  on  density-dependent  models  (Cushing),  it  should 
provide  the  tools  needed  to  begin  the  analysis  of  matrix  models 
and  to  understand  the  relevant  literature. 

Remember  that  the  objective  of  a  matrix  model  (indeed,  any 
structured-population  model)  is  to  describe  the  individual-level 
processes  of  development,  growth,  aging,  survival,  and  reproduc¬ 
tion  to  their  population-level  consequences.  How  much  detail  is 
required  depends  on  the  goals  of  the  model.  To  characterize  popu¬ 
lation  response  to  a  current  environment,  a  linear,  time-invariant 
model  provides  a  powerful  tool.  To  characterize  the  effects  of  popu¬ 
lation  and  environmental  variability,  the  stochastic  tools  in  Chap¬ 
ter  3  are  necessary.  To  include  feedback  between  the  state  of  the 
population  and  the  environmental  conditions  faced  by  the  indi- 
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viduals,  incorporate  density-dependent  nonlinearities  and  use  the 
methods  presented  here  and  in  Chapter  6.  Regardless  of  the  level 
of  detail  in  the  model,  the  analysis  usually  considers  long-term 
dynamics,  problems  of  ergodicity,  and  a  perturbation  analysis. 
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